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ABSTRACT 

We study the effects of anisotropic thermal conduction on magnetized spherical 
accretion flows using global axisymmetric MHD simulations. In low coUisionality plas- 
mas, the Bondi spherical accretion solution is unstable to the magnctothermal insta- 
bility (MTI). The MTI grows rapidly at large radii where the inflow is subsonic. For a 
weak initial field, the MTI saturates by creating a primarily radial magnetic field, i.e., 
by aligning the field lines with the background temperature gradient. The saturation 
is quasilinear in the sense that the magnetic field is amplified by a factor of ^ 10 — 30 
independent of the initial field strength (for weak fields). In the saturated state, the 
conductive heat flux is much larger than the convective heat flux, and is comparable 
to the field- free (Spitzer) value (since the field lines are largely radial). The MTI by 
itself does not appreciably change the accretion rate M relative to the Bondi rate Mb ■ 
However, the radial field lines created by the MTI are amplified by flux freezing as 
the plasma flows in to small radii. Oppositely directed field lines are brought together 
by the converging inflow, leading to significant resistive heating. When the magnetic 
energy density is comparable to the gravitational potential energy density, the plasma 
is heated to roughly the virial temperature; the mean inflow is highly subsonic; most 
of the energy released by accretion is transported to large radii by thermal conduction; 
and the accretion rate M ^ Mb ■ The predominantly radial magnetic fleld created by 
the MTI at large radii in spherical accretion flows may account for the stable Faraday 
rotation measure towards Sgr A* in the Galactic Center. 

Key words: convection, conduction - magnetic fields, MHD - methods: numerical 



1 INTRODUCTION ing local simulations by IParrish fc Stone I l|2005l . 120071 ') and 

„ , r , iParrish fc Quataert I (12003), respectively. 

In magnetized plasmas with a coUisional mean tree path 

much greater than the Larmor radius, there is appreciable In this paper, we study the effects of anisotropic thermal 

heat and momentum transport along magnetic field linos but conduction and the MTI on the global dynamics of spherical 

negligib l e tran sport across the field lines (Braginskii 1965 1. accretion flows using axisymmetric magnetohydrodynamic 

iBalbus I l|2000l 'l found that anisotropic thermal conduction (MHD) simulations. Thermal conduction is particularly im- 

fundamentally alters the convective stability of dilute strat- portant for hot low density plasmas in which the electron 

ified plasmas. For a highly coUisional fluid, convection sets in mean free path is an appreciable fraction of the size of the 

when the entropy increases in the direction of gravity. A low system, as in, e.g., the accretion flow onto Sgr A* in our 

coUisionality plasma, however, is unstable to convective-like Galactic Center and accretion flows onto mass ive black holes 

motions regardless of the sign of the temperature gradient in elliptical galaxies (JTanaka fc Menou 1120061 ). 

(IBalbus II2OO0I : iQuataert 1120081 ). For a temperature increas- mi 1 i- i- r j_ 1 ■ ■ ^ ■ 

^ -^ inere are several motivations tor studying anisotropic 

ing in the direction of gravity, horizontal flelds are the most 1 x- ■ 1 • 1 ^- n n- ^ ^i nj-mr • 

° '=• •" conduction m spherical accretion ttows. lirst, the Mil is 

unstable (the magnetothermal instability; MTI), while when i- j_ j ^ 1 ^ ^ 1 j-- ■ 1 ■ 1 

^ , ' predicted to be present at large radii m spherical accretion 

the temperature decreases m the direction of gravity, vertical „ „ ,, ^»^;2-j_i ■j.i-i 

° •" flows, tor r ^ tb (where rs — GM/c; is the gravitational 
fields are the most unstable (the heat flux driven buoyancy , i- n rj_i ^ 1 1 • ^ 1 ■ j_i • j_i 

, , TT-i-,T\ n Ai r o TT-r,r , , , , Sphere ot luflueuce ot the ccutral oDject aud Cj IS the isotfaer- 

mstabihty; HBI). The MTI & HBI have been studied us- , , 1 r j_i 1 ■ ^ \ a^ j_i i-- ^i 

■' ' mal sound speed ot the ambient gas). At these radii, the 

inflow is subs onic and there is sufficient time for the insta- 
bility to grow (jBalbus 1120001 ) . At smaller radii {r <!^rB), the 

* E-mail: psharma@astro.berkeley.edu MTI is unlikely to grow significantly because the inflow is 
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supersonic (at least according to the classic iBondi Ill952l so- 
lution). However, if there is resistive dissipation of the mag- 
netic energy at small radii, conduction may nonetheless be 
dynamically important because it provides an efficient way 
to transport energy to larger radii, w hich can suppress the 
accre tion rate onto the central object (jJohnson fc Quataert I 
[2007]) • As we discuss in more detail later in the paper, reduc- 
tion of accretion rate does not specifically require anisotropic 
heat transport, but is in fact analogous to the suppression of 
accretion by magnetic dissipation seen in the res istive MHD 
simulations of llgumenshchev fc Naravan I (|2002l ). 

A direct observational motivation for studying the MTI 
in spherical accretion flows is provided by mm observa- 
tions of Sgr A*. These find a rotation measure (RM) of 
~ —6 X 10^ rad m~^ that has been roughly constant since 
it was first measured ov e r ^ 8 years ago ([Aitken et al. I 
I2OOOI : iBower et all l2005l : iMarrone et al. 1 120071 '). The sta- 
bility of the observed RM is difficult to explain if the 
RM is produced by the accretio n flow at small radii 
l|Sharma. Quataert. fc Stone Il2007^ . In local simulations of 
the MTI, the field is amplified by factors of ~ 10 — 
30 and becomes aligne d with the temperature gradient 
l|Parrish fc Stone I [JOOtI ). Thus the MTI at large radii in 
spherical accretion flows is expected to produce a primar- 
ily radial magnetic field that could help explain the stable 
meas ured RM towards Sgr A* (a possib ility we speculated 
on in lSharma. Quataert. fc Stone Il2007^ . 

The remainder of this paper is organized as follows. We 
describe our numerical set up in §2. Because it is not compu- 
tationally feasible to directly simulate both the small radii 
close to the central object where most of the gravitationally 
potential energy is released and the radii ~ rs where the 
MTI develops, we carry out two difi'erent sets of simulations 
in this paper, covering different radial scales in the accretion 
flow. In §3 we present the results of simulations at r ~ rs 
that study the development of the MTI for weak initial mag- 
netic fields (see Table [TJ . In §4 we describe simulations at 
r <^ rB that study the effect of thermal conduction and 
magnetic dissipation on spherical accretion for reasonably 
strong initial magnetic fields (see Table [2| . We then con- 
clude and discuss the astrophysical implications of our re- 
sults (§5). Throughout this paper we consider the simplified 
problem of purely spherical accretion; in future studies, we 
will examine the effects of anisotropic thermal conduction 
and the MTI in rotating accretion flows. 



2 NUMERICAL SIMULATIONS 

We use th e axisymmetric ZEUS-2D MHD code in spherical 
geometry (|Stone fc Norman |[l992al lbh. We solve the MHD 
equations with anisotropic thermal conduction: 



|+V.(pV) = 0, 



(1) 



gX + v.vv^_v(^ + n)+J^_vc, (2) 

at p p 

dB 

'dt 

de 

dt 

Q = -xb(b ■ V)r, 



= V X (Vx B-7]J), 
+ V ■ (eV) = -{p + n)V ■ V - V ■ Q + 7)J^ 



(3) 

(4) 
(5) 



where p is the mass density, V is the fluid velocity, p is 
the pressure, H is an ex plicit artificial viscosity (e.g., see 
ISt one fc Norman |[l992al '). J = c(V x B)/47r is the current 
density, $ = —GM/{r — Vg) is the gravitational potential 
due to a central mass M (we use th e pseudo-Newtonian po- 
tential of lPaczvnski fc Wuta IfigSOD . Vg = 2GM/c^, B is the 
magnetic field strength, r] is the explicit resistivity used in 
some of our simulations, e = p/(7 — 1) is the internal en- 
ergy density (7 is the adiabatic index), Q is the heat flux 
along field lines, x is the thermal diffusivity, b is the unit 
vector along B, and T is the temperature. For convenience 
we often use k = xT/p which has the dimensions of a diffu- 
sion coefficient (cm^ s~^)- Anisotropic thermal conduction 
is a dded in an operator sp lit fashion with subcycling (e.g., 
see IParrish fc Stone 1 12005| ) and is implemented using the 
method based on limiters that guarantees that heat flows 
from hot to cold regions (see Sharma & Hammctt 2007). 

We initialize our simulations using the hydrodynamic 
Bondi solution with a constant vertical magnetic field B^. 
Note that the initial density and temperature gradients are 
in the radial direction; thus, the initial field is not ev- 
erywhere perpendicular to the initial temperature gradi- 
ent, unlike in previous local studies of the MTI. We con- 
sider both weak and strong fields to test the effect of the 
strength of the field on the fiow dynamics. The field strength 
is specified by the dimensionless number Bo, defined by 
Bo = Bz/{8TTGMpout/routy^^, where rout is the outer ra- 
dius of the computational domain and pout = p{rout)- This 



is equivalent to /3(rout) = -Bo [Cs/GM/r 



where (3 is 



the ratio of gas pressure to magnetic pressure. 

We use an adiabatic index 7 = 1.5 so that the initial 
solution has a sonic point that lies within the computational 
domain for our simulations at r ^ r_g (Table [1} [j For our 
simulations at large radii, we use the Bondi radius, r_g = 
GAI/c^, to normalize length scales in the problem, while 
for the simulations at small radii, we use the gravitational 
radius of the central object, rg — 2GM/(? . We express time 
in units of inverse rotation frequency at the outer radius of 
the domain, Q.-^, = [rl^JGHf^''. 

The simulations at ~ tb described in §3 (Table [l| go 
from rin = 0.04rs = 2 x lO^rg to rout = 16rs = 8 x lOVg. 
For these simulations, we keep the temperature at the outer 
radial boundary fixed at its initial value of c?i,/(? — 10"'', 
i.e.. Too i; 10'' K. The simulations described in §4 (Table[2ll 
go from 2rg to 256rg . The temperature is fixed at the virial 
temperature kT^ir/rnp = GMrUp/r at the outer boundary 
(except for one simulation with a very weak initial magnetic 
field [S2]; in this simulation there is no significant change 
in the fiow properties relative to the hydrodynamic Bondi 
solution, so fixing the temperature to be that of the ini- 
tial Bondi solution is more physical). For the simulations 
at ^ rs (§3), fixing the temperature at the outer bound- 
ary is reasonable because of the large thermal inertia of the 
plasma at large radii. For the simulations from 2rg to 256r5 



^ For 7 = 5/3 the sonic point lies very close to the origin. It is, 
however, numerically desirable to have the sonic point within the 
computational domain so that the inner boundary conditions do 
not influence the solution. We have verified that we get similar 
results for a run with 7 = 5/3, but with rest of the parameters 
same as the fiducial run Rl. 
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Table 1. Simulation parameters for runs with r ^ tb (§3) 



Name 


Resolution 


4 


Bl 




Rl" 


60 X 44 


0.2 


4.5 X 10" 


-4 


R2 


60 X 44 


0.2 


4.5 X 10" 


-5 


R3 


60 X 44 


0.2 


1.34 




Dl 


120 X 88 


0.2 


4.5 X 10- 


-4 


D2 


120 X 88 


0.2 


4.5 X 10" 


-5 


Ql 


240 X 176 


0.2 


4.5 X 10" 


-4 


MHDl 


60 X 44 





4.5 X 10- 


-4 


MHD2 


60 X 44 





1.34 





Hn = 0.04rB, rout = 16rfl 

No explicit resistivity is used for these simulations. 

'^ This is our fiducial MTI simulation. 

t The conductivity is defined by k{t) = ac{GMr)^^^ (see il2.1| l. 

t Bo is the dimensionless magnetic field strength, given by Bo = 

B^/(87rGMpout/r-out)^/2. 



Table 2. Simulation parameters for runs with r ^ tb (§4) 



Name 


Resolution 


4 


Bl 


Mb 


SI 


60 X 44 


0.2 


0.32 


0.01 


S2 


60 X 44 


0.2 


10-" 


0.52 


S3 


60 X 44 


2 


0.32 


0.009 


S4 


60 X 44 


20 


0.32 


0.039 


MHD3 


60 X 44 





0.32 


0.028 



rin = 2rg, rout = 256rg 

An explicit resistivity is used for these simulations (eq. [6]). 

t The conductivity is defined by K(r) = ac{GMr)^''^ (see ^2.11 1. 

f Bo is the dimensionless magnetic field strength, given by 

So = -B^/(87rGMpout/rout)i/2. 

* time averaged accretion rate through rin. 



described in §4, however, the choice of the outer tempera- 
ture boundary condition is more subtle and has a nontrivial 
effect on the solution. Because we initialize the flow with 
a 7 = 1.5 < 5/3 solution, the initial temperature is sig- 
niflcantly less than the virial temperature. Fixing Tout to 
its initial value is inappropriate because dissipation leads to 
a temperature ~ Tvir throughout most of the domain and 
thus setting Tout equal to its initial value leads to an unphys- 
ically large temperature gradient at rout- Setting dT/dr — 
is unphysical because this sets the heat flux to 0. An al- 
ternative boundary condition that we considered, setting a 
constant heat flux across the outer boundary (corresponding 
to dT/dr ~ constant), led to a roughly constant tempera- 
ture throughout the domain with Tout — Tn ^ Tvir after 
several fiouti this drove all the mass out of the simulation 
domain. We thus settled on the physically motivated outer 
boundary condition of Tout = Tvir- Factor of few variations 
in the exact value of Tout did not significantly change the 
results of the simulations. In both sets of simulations, in- 
flow/outflow boundary conditions are applied at the radial 
boundaries. 

The radial grid is logarithmically spaced (grid spac- 
ing dr oc r) so that the number of grid points from rin to 
(j'inrout)^ is roughly the same as the number of grid points 



from (nn^out)^ to rout- A logarithmic grid gives reasonable 
resolution at all radii; e.g., for our simulations at small radii 



uniform polar grid extends from ^ = to ^ = vr. Reflective 
boundary conditions are applied at the poles. Our standard 
resolution is 60 x 44, but we have also carried out simulations 
at higher resolution to assess convergence. Our calculations 
use a version of ZEUS that is not parallel, so that even these 
modest 2D calculations are time intensive. 

We do not use an explicit resistivity for the r "^ tb 
simulations in §3 (Table [T|. These simulations are primarily 
aimed at studying the MTI for weak fields in which case 
the magnetic energy is so small that using an explicit re- 
sistivity is not required to correctly capture the dynamics 
(in addition, we have verified that using the explicit re- 
sistivity below does not change the results). However, for 
r <^rB (§4; see Table [2)), the gravitational potential energy 
of the inflowing plasma is converted into magnetic energy 
by flux freezing. The infiow brings together oppositely di- 
rected field lines, which leads to magnetic dissipation and 
plasma heating. This energy is then transported to large 
radii by thermal conduction. To capture this dynamics, it 
is important to conserve energy as accurately as possible. 
Although we are not using a conservative code, we can cap- 
ture a reasonable fraction of the dissipate d magnetic energy 
using an explicit resistivity of the for m l|Stone fc Pringle I 



using an explicit resistivity ot the tor r 
l200ll : llgumenshchev fc Naravan II2002I ) 



rj = 770 dr' 



IV X Bl 



(6) 



with rjo = 0.15 (with this choice for 77, we find that total 
energy is conserved to better than 30 % in our simulations 
at r <C vb)- The resistive terms in the induction and in- 
ternal energy equations are inclu ded using the method of 
[Fleming. Stone, fc Hawlev I (|2000| ). 

Although our simulations are two-dimensional, the tur- 
bulent magnetic and velocity fields do not decay at late 
times, because the anti-dynamo theorem does not apply 
when a net magnetic fiux threads the simulations domain. In 
this case the source curre nt for the magnetic field i s outside 
of the simulation domain (jBalbus fc Hawlev Ill998l ). 

2.1 Thermal Conductivity 



In a coUisional plasma, the thermal conductivity is given 
by K ~ 4.3 X lO^T^/V"' cm^ s"^ l|Spitzer I [l963 ) . How- 
ever, the inner regions of hot accretion flows are, in 
many cases, collisionless with the electron mean free 
path due to Coulomb collisions larger than the radius 
(e.g., iTanaka fc Menou I l2006l : ISharma et all |2007| ): the 
coUisional result is thus inapplicable. In a collisionless 
plasma, electrons can in principle carry a free stream- 
ing heat fiux as large as Q ~ n kTve ~ 40/9tif (v j. 



jkTl 



11/2 



is the ther mal spcod JCowie fc McKee Ill977l : 



ISnvder. Hammett. fc Dor land 19971). The diffusion coeffi 
cient in this case is k ~ Ver oc r^'^, assuming a virial tem- 
perature proportional to r~^ . Motivated by the application 
to collisionless systems such as Sgr A*, we take n oc r^'^ as 
the thermal conductivity for our simulations and define the 
dimensionless conductivity ac by 

K = ac{GMrf^^. 



(7) 

Although the scaling k ~ VeV oc r^ " is well-motivated, the 
precise value of Oc is difficult to calculate. It depends on the 
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electron mean free path in the coUisionless hmit, which is de- 
termined by the rate of pitch angle scattering by small-scale 
(high frequ ency) turbule nt fluctuations. Motivated by the 
results of S harma et al. I (|2006 ) and Sharma et al. I (|2007l ') 
on pitch angle scattering in coUisionless accretion flows, we 
take Oc ~ 0.2 - 2, i.e., k = (0.2 - 2)(G'Mr)i/^ as our fidu- 
cial diffusion coefficient. This value is much smaller than the 
maximum free-streaming thermal conductivity that could 
be obtained if the electrons are virial. Our choice for «;(r) 
is such that the large-scale thermal conduction timescale 
(tcond = r"^ / 1"^) is comparable to the free-fall timescale 
(tft = [2GA//r^]~^/^). Thermal conduction can thus lead 
to a significant modification of the thermal structure of the 
plasma (as we shall see). Since both icond and tg scale as 
r^", this is true at all radii. In addition, with icond ~ iff on 
large scales, small-scale perturbations are effectively isother- 
mal and the MTI grows at its maximal rate (Balbus 2000|). 
For the MTI simulations at ~ rs , we carried out simulations 
with Uc factors of few smaller and larger than our fiducial 
choice of Qc = 0.2, and found results essentially identical to 
those described here. For the simulations at r ~ 2 — 200rg , 
resistive heating and thermal conduction significantly mod- 
ify the dynamical properties of the accretion flow; in this 
case, we show solutions for a variety of values of Qc to ex- 
plicitly show the dependence on the conductivity (which is 
still relatively modest). 

As we discuss in §3, it turns out that a conductivity 
K oc r^" has the property that V • Q ~ for the initial 
Bondi solution near r ^ rs- Since the initial condition is 
already approximately a steady state solution even in the 
presence of conduction (modulo the MTI), including con- 
ductivity does not affect the structure of the plasma. We 
thus carried out simulations at ~ rs with other choices for 
K{r) (namely, k cx r^'^ and k cx r~^'^). Although we do 
not show detailed results for these simulations, in all cases 
we found that the MTI rapidly creates radial magnetic field 
lines and that conduction along the radial magnetic field 
lines forces the density and temperature structure of the 
plasma to readjust in order to satisfy V • Q ~ 0. This can 
lead to modest modifications in the density and temperature 
profiles of the fiow, but in all cases the accretion rate M was 
within a factor of a few of the Bondi rate. 



3 SIMULATIONS AT i? ~ i?s 

Figure[T]shows several different timescales characterizing the 
initial Bondi solution. For r > vb the infall timescale (tin) is 
much longer than the timescale for the MTI to grow (iMTi) 
and the instability grows exponentially. The MTI can thus 
saturate before the plasma flows in. At smaller radii (r < 
rs), however, where the inflow time is comparable to the 
MTI growth time, the growth is algebraic and the MTI is less 
likely to have a significant effect on the dynamics. Instead, 
at r ^ r_B, the magnetic field will be primarily amplified by 
fiux freezing. Here we discuss our simulations of the MTI 
at r ~ Tfl. In the next section we present the results of 
simulations at r <^rB- 



a 



1000 FT 

100 r 

10 r 

1 r 

0.1 r 

0.01 r 

0.001 ^ 



t 




0.0001 



Figure 1. Various timescales (in units of C~^(. = [GM/r^^^]"^''^) 
as a function of radius for the initial Bondi solution: the MTI 
growth time tMTl = 7mti (*^^ maximum MTI growth rate is 
given by 'yl^rj,^ = [l/p][dp/dr][dlnT/dr]), tin = r/Vr (the local 
infall timescale), t^ond = r-^/^, (the large scale thermal diffusion 
timescale for ac = 0.2), and fff = {r^/2GMy/^ (the local free-fall 
timescale). The infall timescale is much longer than the timescale 
for the growth of the MTI for r > rs] thus the MTI rapidly 
modifies the magnetic field structure. 



3.1 The Fiducial MTI Simulation 

Our fiducial simulation (labeled "Rl") uses K(r) = 
0.2(GAfr)^" and has an initial magnetic field strength of 
Bo = 4.5 X lO""*, corresponding to /? ~ 10* at rout. This 
simulation describes the evolution of the MTI for weak mag- 
netic fields. We choose the lower resolution simulation Rl as 
our fiducial case because it finishes relatively quickly (in < 
a week); this enables us to readily study the effects of dif- 
ferent boundary conditions, different choices for the thermal 
conductivity, etc. We also present higher resolution calcula- 
tions (e.g., Dl & Ql, which take several months to finish) 
that are similar to Rl, although with a modest increase in 
the field amplification by the MTI (see Fig. [S]). 

Figure [2] shows the r and 9 components of the mag- 
netic energy density as a function of time for the fiducial 
simulation and for an MHD simulation without conduction, 
but with the same initial conditions. The energy is aver- 
aged from r — rout/10 to rout, to isolate the radii where 
tin 3> tiviTi (see Fig.[l| and the instability is unaffected by in- 
flow. The top panel in Figure[2]shows the early linear growth 
of the MTI, while the bottom panel shows the longer time 
evolution. Because we are averaging over a range of radii, 
there is no unique growth rate for the MTI, but the maximal 
MTI growth rate evaluated at r ~ rout/3 provides a reason- 
able fit to the early time growth (long-dashed line). After 
a few free-fall times (several MTI growth times; see Fig. 
[l}, the magnetic field ceases to grow exponentially and the 
growth becomes algebraic. This represents the saturation of 
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10 e 




0.001 



^o.tt 




r, R] 

0, Rl - 

r, MHDl- 



Figure 2. Change in volume averaged magnetic energy (radial 
and 6 components normalized to the initial magnetic energy) in 
a shell from rout/10 to rout as a function of time. The top panel 
shows the early-time evolution of SB^ (solid line) for the fidu- 
cial MTI simulation (Rl) compared to SB^ for an MHD simula- 
tion with the same initial conditions (MHDl; short-dashed line). 
The long-dashed line indicates the growth rate of the MTI at 
r ~ rout/3. The bottom panel shows the same field components 
over longer timescales, along with SBg (dotted line) for the fidu- 
cial run. SB^ g = ({B'^g)v — B^ g q)' where Qy denotes volume 
average and Bq is the initial field strength. 



the MTI; the later growth of the field is due to flux freezing. 
This can be seen explicitly by noting that the MHD simula- 
tion (short-dashed line) shows the same algebraic growth in 
the magnetic field with time, only at a reduced magnitude 
of B because there is no early MTI growth. Note that the 
net increase in magnetic energy because of the MTI is quite 
modest: /3 3> 1 even in the saturated state. 

Figure[3]compares the plasma (3 at the end (flouti — 60) 
of the fiducial run and the corresponding MHD simulation. 
The magnetic field unit vectors are also shown. The plasma 
/3 with conduction is a factor of ~ 100 smaller than without 
conduction (as is also seen in Fig. [2]), indicating that the 
MTI alone amplifies the field by a factor of ~ 10. Since the 
simulations do not run for even one infall timescale at rout, 
the magnetic field lines at large radii do not change signifi- 
cantly in the MHD simulation and remain primarily vertical. 
By contrast, the magnetic field is completely restructured in 
the presence of anisotropic conduction because of the MTI. 
The magnetic field at the end of the simulation is primarily 
radial. This can also be seen in Figure ID which shows the 
volume averaged angle of the magnetic field with respect to 
the radial direction as a function of time for several runs. 
In the hnear stage of the MTI, the radial field is amplified 
exponentially; after saturation, the field direction, although 
turbulent, continues to become radial because of fiux freez- 
ing (as also occurs in the pure MHD simulation shown in 
Fig. |4|, although this occurs on the slower inflow timescale. 




I 






Figure 3. Logj^Q of plasma /9 (= Sirp/B^) at the end (f2outi ~ 60) 
of the fiducial MTI run Rl (left), and at the end of the MHD sim- 
ulation MHDl (right). Also shown are arrows indicating the local 
magnetic field direction. The magnetic field is turbulent but pri- 
marily radial in the presence of anisotropic conduction. In MHD, 
the field remains primarily vertical at large radii, because fiux 
freezing only reorients the field on the (slow) inflow timescale. 
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Figure 4. Volume averaged angle (in degrees) of the magnetic 
field lines with respect to the radial direction as a function of time 
for difl'erent runs (the average is from rout /ID to rout, as in Fig. 

Ell. 



Our conclusion that the MTI generates primarily radial fleld 
lines in the saturat ed state is fully consistent with the local 
MTI simulations of iParrish fc Stone I (|2005l . l2007l ). 

To understand the saturation of the MTI, it is useful to 
consider the linear dynamics of the instability. In the WKB 
limit, and assuming rapid conduction and weak magnetic 
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fields, the Unear growth rate is given by (JQuataert Il2008f ) 



•y Si 



I dp f dlnT 
p dr 



dr 



,^ _ r.i'Z-.k'g 2bgbrk0kr 



100 



fc2 



(8) 



where kr and kg are the r and 9 wavenumbers, and br and bg 
are the r and 9 magnetic field unit vectors. Our initial field is 
vertical with a 9 dependent br and bg. In the saturated state, 
where 6r — 1 and bg <^1, the fastest growth rate of the MTI 
is 7 ~ \\1/ p\[dp/dr\[d\nT/dr\\^''^bg, a factor of ~ bg smaller 
than the growth rate in the initial configuration. Growth 
occurs only for modes satisfying kg/kr < 2bg, i.e., for small 
radial wavelengths. The significant decrease in the growth- 
rate of the MTI for nearly radial fields qualitatively explains 
why magnetic field reorientation leads to saturation of the 
instability. In addition, even if there is growth at late times, 
it is on small radial scales and so is unlikely to modify the 
large-scale field structure or magnetic energy. The MTI can 
also saturate by making dlnT/dr « 0, i.e., by making the 
plasma isothermal. Whether the plasma becomes isothermal 
or not depends on K{r) (and on the boundary conditions), 
which infiuences how the density and temperature adjust 
to make V ■ Q ~ in steady state. Figure [5] shows that 
for our fiducial simulation, there is very little change in the 
temperature and density profiles during the simulation. For 
our choice of k oc r^' ^, it turns out that the initial Bondi fiow 
nearly satisfies xr^dT/dr — constant near ~ tb- Thus when 
the MTI reorients the field lines to be primarily radial, little 
evolution in p{r) and T{r) is required to make V ■ Q ~ 0. As 
noted in i]2.1l we carried out simulations with different fi:(r) 
in order to assess its effects on the flow dynamics (although 
we believe that the fiducial choice considered here is the most 
physical for collisionless systems). In all cases, we found that 
the MTI reorients the field lines to be primarily radial and 
that p(r) and T{r) adjust so that V • Q ~ 0. 

Figure [6] shows the time and angle averaged mass ac- 
cretion rate in the saturated state as a function of radius. 
In MHD the velocity is always radial and inwards and the 
mass accretion rate is equal to the Bondi rate (for this /3 3> 1 
simulation) . The solid line in Figure [S] shows the net mass 
accretion rate for the fiducial run, M = 2-Kr^ J pVr sin 9d9; 
also shown are the mass fiuxes of material with V,- > and 
Vr < 0. Although there are MTI-driven convective motions, 
they are modest and the net accretion rate is only slightly 
smaller than the Bondi accretion rate. This is partially a 
consequence of the fact that the density and temperature 
profiles change very little in the course of the simulation 
(Fig. [5]). However, even for different choices of Oc and «:(r), 
we still found M ~ Mb for MTI simulations with initially 
weak fields. The MTI by itself thus does not significantly 
change the accretion rate onto the central object. This is in 
contrast to the very strong suppression of M found in the 
next section by the combined action of magnetic dissipation, 
plasma heating, and thermal conduction at r ^ rs. 

Figure [7] shows time and angle averaged luminosities 
(47rr^x energy fiuxes) for energy transport by conduction 
and MTI-initiated convection. The luminosities are normal- 
ized to the field free ("Spitzer") conductive luminosity at 
rout- The radial conductive luminosity (Airr'^Qc) is roughly 
constant as a function of radius for r ^ O.Irs, as would be 
expected in steady state. The conductive heat flux is similar 
to the local field-free value at all radii because the field lines 
are primarily radial. The non-conductive energy transport 
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Figure 5. Density (short-dashed line) and temperature (long- 
dashed line) as a function of radius for the initial Bondi solu- 
tion; time and angle averaged density (solid line) and tempera- 
ture (dotted line) at late times as a function of radius for the 
fiducial MTI simulation Rl. Density and temperature are nor- 
malized to their values at rout. There is very little change in the 
radial density and temperature profiles during the simulation. 
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Figure 6. Time and angle averaged mass accretion rate {M) in 
the fiducial MTI simulation Rl at late times, as a function of 
radius (normalized to the Bondi accretion rate, Mb)'- uet mass 
accretion rate (solid line), mass inflow rate (dotted line; with Vr < 
0), and mass outfiow rate (short-dashed line; with Vr > 0). The 
mass accretion rate in MHD is shown by the long-dashed line. 
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Figure 7. Time and angle averaged conductive (solid line; 
Qc) and turbulent (dotted line; Qt) luminosities (Anr-^ X energy 
fluxes) as a function of radius for the fiducial MTI simulation Rl; 
see §3.1 for the definition of Qt- The luminosities arc normalized 
to the field free conduction luminosity at rout- 



Figure 8. Volume averaged magnetic energy (normalized to the 
initial magnetic energy) in a shell from rout/ 10 to rout as a func- 
tion of time: SB^/B^ (solid line), SBJ/B^ (dotted line) for the 
fiducial MTI simulation Rl; SB^/B§ (short-dashed line), SB'^/B^ 
(long-dashed line) for R,2 (which has an initial field strength So 
10 times smaller than Rl); SB^, g = {{B^ g)Y — B^ g „), where Qv 
denotes volume average and Bq is the initial field strength. 



has contributions from the kinetic energy flux {pV'^Vr/2), 
the enthalpy flux (7/(7 — l)pVr), and the Poynting flux. 
The latter is negligible for the weak magnetic field values 
considered here. Both the kinetic energy flux and the en- 
thalpy flux can be split up into terms proportional to the 
mass accretion rate (the advective flux of energy Qa due to 
the bulk motion of the fluid) and terms which represent the 
turbulent/convective transport of energy (Qt). Mathemat- 
ically, this can be seen by decomposing the fluid variables 
into mean and fluctuating components; e.g., the radial ve- 
locity is taken to be Vr = (Vr) + SVr, where () represents 
an appropriate average and S represents deviation from the 
average. For our purposes, the average is taken over 6 at 
each time at a given r. With this decomposition, the aver- 
age kinetic energy flux is given by 



-{pV'Vr)=QK,a + QK,t, 



with 



QK,a = ^{pVr)[{Vrf + {Vof + {SK') + {SVi)], (9) 



QK,t = -{5{pVr){5V; + SVe')) + {S{pVr)SVr){Vr) 

+ {5{pVr)5Ve){Ve), 



(10) 



where Qx,a is the advective flux of kinetic energy and QK,t 
is the turbulent/convective kinetic energy flux. Similarly, the 
average enthalpy flux can be written as 



7 



7-1 



{pVr) = Qe,a +Qe.t, 



with 



^e,a 



^ {pVr){cl), 



7-1 



= ^ liVr) {5p5cl) + (p) {SVrScl) 

7- 1 
-I- {Sp5Vr&cl)], 



(11) 



(12) 



where c^ = p/ p ~ kT/rUp. Figure [7] shows the total tur- 
bulent energy flux in the MTI simulations at large radii, 
Qt = Qe,t + QK,t\ the turbulent energy transport - which is 
dominated by the {SVrSc?^) term in equation (|12|) - is much 
smaller than the conductive energy transport at all radii. 



3.2 Effect of Bo 

If the MTI saturates by making the magnetic fleld lines 
primarily radial, as we showed in the previous section, it 
is reasonable to expect that the saturation magnetic fleld 
strength will scale with the initial magnetic fleld strength 
(since a smaller radial field is required to make an initially 
weaker vertical magnetic field radial). To test this hypothe- 
sis we carried out the simulation R2 with a weaker vertical 
magnetic field {Bq = 4.5 x 10~^), corresponding to /3 ~ 10^" 
at Tout. All other parameters are same as the fiducial run. 

Figure [S] shows the results from run R2 compared to 
those of Rl. The most interesting result is that even nonlin- 
early SB^ g/Bo for Rl and R2 are similar. That is, the late 
time magnetic energy is roughly proportional to the initial 
magnetic energy in the domain. This indicates that the sat- 
uration of the MTI is quasilinear: for weak magnetic fields. 
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the MTI saturates when the field fines become radial, wfiicfi 
requires a SB ex Bo. 

We also carried out simulations with mucfi stronger ini- 
tial magnetic fields, Bo = 1-3, corresponding to /3 ~ 10 at 
''out (run R3 with conduction and run MHD2 without con- 
duction). The magnetic field is sufficiently strong that the 
MTI is suppressed by magnetic tension. Indeed, we find that 
there are only small differences between the simulations with 
and without conduction. In both cases the magnetic field 
at small radii is amplified by the converging inflow. Mag- 
netic tension prevents plasma infall in the equatorial region; 
mass accretion is dominated by the polar regions. There 
is a strong equator-pole anisotropy because of the strong 
magnetic field. The polar regions are cold, low density, and 
magnetically dominated (/3 ^ 1), with large radial infall ve- 
locities. We find a modest reduction in the accretion rate 
for these strong-field simulations at ~ rs, with the mass 
accretion rate through rin being ~ 25% of the Bondi rate. 
Finally, we find that including an explicit resistivity does not 
have a significant effect on the flow dynamics at large radii 
even for energetically signiflcant magnetic fields. This is in 
contrast to the simulations described in §4 in which there is 
significant magnetic dissipation when /3 ~ 1 at r <C r_B . The 
difference is probably that the supersonic inflow at r <C rs is 
much more effective than the slow inflow velocities at ~ tb 
at forcing reconnection/dissipation of fleld lines. 



10'' 



3.3 Effect of Resolution 

To assess the effect of resolution we carried out several higher 
resolution runs: run Dl, which is analogous to the flducial 
simulation Rl but with twice the resolution, run D2 which 
is analogous to R2 but with twice the resolution, and run Ql 
which has four times the resolution of Rl. While the overall 
dynamics is very similar to the lower resolution runs, the 
magnetic energy is a factor of few larger in the higher res- 
olution simulations. Figure [9] compares the flducial run Rl 
with Dl and Ql, and shows that both 5B^ and 5B^ are 
larger in higher resolution runs than they are in Rl, by a 
factor of ~ 5. The same is true for D2 compared with R2. 
In particular. Figure [5] shows that the phase of linear ex- 
ponential magnetic fleld ampliflcation occurs for a longer 
period of time in the higher resolution simulations. Since 
the late-time magnetic field growth is simply fiux freezing 
of the early-time field, this leads to an enhanced field at all 
times. The fastest growing MTI modes for approximately 
radial magnetic fields are those with 2fcr/fce >bg^, which is 
a resolution dependent criterion. Higher resolution simula- 
tions can thus resolve growth for longer times and to larger 
radial fields, i.e., to smaller bg. 

Although the fastest growth of the MTI occurs at small 
scales, nonlinearly one expects the energy contained in small 
scales to be subdominant for sufficiently high resolution sim- 
ulations. Indeed, Figure [9] shows that both the radial and 
polar magnetic energies are similar for the higher resolution 
runs Dl and Ql. This indicates that nonlinear saturation 
of the MTI is reasonably independent of resolution for suf- 
ficiently high resolution simulations (e.g., runs Dl and Ql). 
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Figure 9. Volume averaged magnetic energy (normalized to the 
initial magnetic energy) in a shell from rout/ 10 to rout as a func- 
tion of time, comparing different resolution runs: 5B^/Bq (solid 
line), 5Bl/Bl (dotted line) for the fiducial run Rl; SBI/BI 
(short-dashed line), SB^/Bq (long-dashed line) for Dl, with 
twice the resolution as Rl; and SB'^/Bq (short dot-dashed line), 
SB^/Bq (long dot-dashed line) for Ql, with four times the res- 
olution of Rl. As before, SB^g = {(B^g)Y — B"^ g p), where Qy 
denotes volume average and Bq is the initial field strength. 



4 SIMULATIONS AT R ^ Rb 

As shown in the previous section, the MTI generates a ra- 
dial fleld at ~ rB in dilute spherical accretion flows. As 
the plasma flows in to yet smaller radii, this fleld will be 
amplifled by flux freezing. It is well known that for inflow 
according to the Bondi solution (Bondi 1952 ), the magnetic 
energy increases more rapidly with decreasing radius than 
the gravitational energy of the gas, an d eventually the flel d 
becomes energetically imp o rtant (e.g., IShvartsman I Il97ll ) . 
Ilgumenshchev fc Naravan I (|2002j showed that when /3 ~ 1 
magnetic dissipation of the amplified field becomes signifi- 
cant and the dynamics of the infiow changes appreciably. Re- 
sistive heating increases the temperature of the plasma and 
may drive thermal convection. The accretion rate onto the 
central object decreases by a few orders of magnitude. In this 
section, we examine the modification to these results caused 
by anisotropic thermal conduction. Our general conclusion 
that mass accretion is significantly suppressed is consistent 
with that of Igumcnshchov & Naravan (2002), but we find 
that the conductive transport of energy is more important 
than the convective transport. 

Our simulations at small radii have rin = 2rg and 
^out = 256rg ; their properties are summarized in Table [51 
The initial conditions are chosen from the same Bondi so- 
lution used in the previous section. Run SI takes Qc = 0.2 
and Bq ~ 0.32. This value of B q corr esponds to /? ~ 1 at 
Tout. Ilgumenshchev fc Naravan I (|2002h have shown that a 
weaker field at rout means that significant dissipation does 
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Figure 10. Mass accretion rate (normalized to the initial Bondi 
rate) at ri^ = 2rg as a function of time for runs SI, S2, S3, 
& MHD3, all listed in Table [2] The mass accretion rate in the 
presence of strong fields and resistive dissipation is reduced sig- 
nificantly relative to the Bondi rate, irrespective of whether con- 
duction is present. Note that the mass accretion rate for run S3 
shows very regular "flares;" these persist even if we vary many of 
the numerical details of the calculation (see text for details) . 
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Figure 11. Total energy generation and loss rates (normalized to 
GMMfl/rin) as a function of time for run S3 with Oc = 2: gravi- 
tational energy released by accretion GMM/ri^ (solid line), vol- 
ume integrated resistive heating rate (H = 2tt ^ drdGr^ sinOr)J^\ 
dotted line), and conductive luminosity at the outer radius of 
the domain (47rr^Qc|out; short-dashed line). Most of the energy 
released by accretion is captured by our explicit resistivity; this 
energy is then transported to large radii by conduction. 



not begin until smaller radii (and later times) , but that oth- 
erwise the results do not depend sensitively on the choice 
of Bq. In addition to the fiducial value of a^ = 0.2, we also 
performed simulations with Oc = 2 and 20; the latter corre- 
sponds roughly to the maximum rate of heat transport by 
free-streaming electrons. We present results for a variety of 
ac but focus on the Qc = 2 simulation (S3) because it most 
clearly demonstrates the relevant physics. 

S2 and MHD3 represent control simulations. S2 is a sim- 
ulation with anisotropic conduction but with a very weak 
field. As expected, it is effectively a laminar hydrodynamic 
Bondi solution (see Fig. ll4|l . with no evidence for the MTI or 
dynamically significant resistive heating. MHD3 is a resistive 
MHD simulation wit h /3 ~ 1 at rput, an alogous to the simu- 
lations performed by llgumenshchev fc Naravan I (|2002h . 



4.1 Results 

Figure [10] shows the mass accretion rate at the inner ra- 
dius rin as a function of time. The weak-field simulation S2 
shows no significant differences relative to the Bondi solu- 
tion; the factor of ~ 2 decrease in M for S2 in Figure [10] 
is due to the fact that thermal conduction leads to small 
changes in the density and temperature profiles of the fiow. 
In addition, there is no indication that the MTI is present, 
unlike in the simulations at large radii discussed in the pre- 
vious section. By contrast, all of the simulations with /3 ~ 1 
at rout - both with and without anisotropic conduction - 
show a significant decrease in M after a few free-fall times. 



with M ~ 0.01 — 0.03 A'Ib at late times. We focus our inter- 
pretation on the simulations with conduction, since our re- 
sults for simu lations without conduction a re sim ilar to those 
presented in llgumenshchev fc Naravan I (|2002i ). Note that 
the simulations with anisotropic conduction show episodic 
"flares" in M and other dynamical quantities that are not 
as prominent in the MHD simulations; these are discussed 
more below. 

Figure[Tl]shows the overall energetics for the simulation 
with cic = 2 (S3). A significant fraction of the gravitational 
potential energy released by the inflowing matter (solid line) 
is converted into magnetic energy by flux freezing. Resistive 
dissipation of oppositely directed field lines leads to plasma 
heating. Although our numerical routine is not conserva- 
tive, we do have an explicit resistivity and Figure [TTI shows 
that the total resistive heating in the computational domain 
(dotted line) is similar to the gravitational potential energy 
released by accretion[j In turn, much of this energy is trans- 
ported out of the computational domain by conduction, as 
shown by the short-dashed line in Figure 1111 which is the 
conductive luminosity at the outer radius rout- The results 
for runs SI & S4 (which have different conductivities Oc) 
look similar. 



^ In addition to explicit resistivity — which is the dominant source 
of heating in our calculations — we also include artificial viscosity 
for shock capturing. Explicit resistivity, artificial viscosity, and 
inherent numerical diffusion all lead to increase in entropy in the 
simulations. 
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Figure 12. Time and angle averaged luminosities (47rr X en- 
ergy fluxes, normalized to GMMb/i"\d) as a function of radius 
for run S3 (qc = 2): conductive luminosity iirr^Qc, energy trans- 
port by fluid motion 47rr'^Qe+x-|-s (which includes thermal, ki- 
netic, and Poynting contributions), total power generated by ac- 



cretion (GMM/[r - Tg]), 
where Qtot =Qc + Qe+K- 



and the total luminosity (Attt Qt. 

B + GMM/i-n-lr - r^]). 



Another indication of the importance of the conductive 
transport of energy is given in Figure 1121 which shows the 
time averaged luminosities carried by conduction (Qc) and 
fluid motion (Qs+k+b) as a function of radius for etc = 2 
(S3). Note that here we are comparing the conductive lu- 
minosity to the total thermal-|-kinetic-|-magnetic transport 
of energy Q,+k+b = pV^Vr/2 + 7PK/(7 - 1) + [B^Vr - 
(B • V)Vr]/47r which includes kinetic, thermal, and Poynt- 
ing contributions; the radial velocity Vr in Qe-i-K+s is the 
total radial velocity and thus Qe+x-i-s includes energy trans- 
port by the mean fluid motion and by the turbulence (unlike 
Qt in §3 and Fig.[7]which is solely the turbulent transport). 
The energy carried by conduction is significantly larger than 
that carried by other mechanisms (in particular, either ad- 
vection or convection) except at the smallest radii, where 
the pseudo-Newtonian potential leads to supersonic inflow. 
In addition, the conductive luminosity is roughly indepen- 
dent of radius, as would be expected in steady state given a 
source of heating at small radii. These results demonstrate 
that the gravitational potential energy released at small radii 
by magnetic dissipation is carried to larger radii by conduc- 
tion. We find similar results for Oc = 0.2 (SI) and ac ~ 20 
(S4), but for larger values of Oc, the conductive energy trans- 
port becomes the dominant energy transport mechanism 
at smaller radii. Only in the limit ac ^ 1 do the simula- 
tions with conduction become sim ilar to the simulations of 
llgumenshchev fc Naravan I l|2002r ). 

Figure [T5] shows how the time and angle averaged den- 
sity and temperature profiles of the accretion flow change in 
the presence of resistive heating and conduction. The initial 



temperature of the flow is much less than the virial tem- 
perature fcTvir = GMnip/r because we initialize a Bondi 
solution with 7 = 1.5 < 5/3. By contrast, both the MHD 
simulation and the simulations with anisotropic conduction 
have time-averaged temperatures of T ~ Tvir because much 
of the accretion energy is thermalized as heat via resistive 
dissipation. The resulting high temperatures lead to a flow 
in approximate hydrostatic equilibrium, with a mean radial 
velocity that is highly subsonic at r > a few rg (see Fig. ll4p . 
Figure [13] also shows that the density of the flow in the inner 
region is much lower in the simulations with heating than 
in the initial solution. This is consistent with accretion rate 
being < Ms (Fig. fTO)) . 

In Figures 112)1141 we have presented results averaging 
over the full range of 0, from 9 — to 9 — n. However, 
the initially vertical magnetic field lines impose an asymme- 
try on the flow so that the flow profiles are not spherically 
symmetric. As seen in Figure \W\ (discussed below), the tem- 
perature is larger in the equatorial region compared to the 
pole, by roughly a factor of 10 at small radii; the same is 
true for the density. The equatorial region is primarily gas 
pressure dominated while the poles are magnetically domi- 
nated. The relative importance of the different energy trans- 
port mechanisms (shown in Fig. I12p also varies with angle; 
in particular, conduction dominates at the poles (and aver- 
aged over all angles) , but is less important in the equatorial 
region where the field lines are less ordered. Overall, the an- 
gular asymmetry we find in these nonrotating simulations is 
not as large as in simulations with s ignificant angular mo- 
mentum (e.g.. lHawlev fc Balbusll2002r ) which show "funnels" 
with radial outflows. In our spherical simulations, the time 
averaged flow is always inwards at all angles (although there 
can be outflow at certain times; see Fig. I16|l . 

Although there are modest quantitative differences be- 
tween the simulations with anisotropic conduction and those 
without (Figs. I10m4p . our results demonstrate that many 
of the properties of magnetized spherical accretion are in 
fact rather similar in the two cases. An important differ- 
ence, however, is that for even modest thermal conductivities 
(qc ~ 0.2; eq. [7]), conduction is the dominant mechanism of 
energy transport and transports a significant fraction of the 
energy released by magnetic dissipation to large radii (Fig. 
llll fc [T2)l . By contrast, in the MHD simulations without con- 
duction, we find that the net energy transport is inwards. 
Although convective motions are present, the convective lu- 
minosity is much smaller than GMM /{r — Vg) at all radii. 

An additional interesting difference between the simu- 
lations with and without conduction is the "flaring" in the 
simulations with anisotropic conduction (Fig. I10lll|l . The 
accretion rate, resistive heating, and conduction luminos- 
ity all show regular rapid increases ( "flares" ) by a factor of 
~ 10 — 100 on a timescale ~ 10 orbital periods at ~ rin (see 
Fig. I15|) . The characteristic timescale between such flares is 
a fraction of an orbital period at rout. The flare timescales 
do not depend sensitively on the value of the thermal con- 
ductivity, as long as a ~ 2. As shown in Fig. 1101 flares are 
most prominent for run S2, the run with the smallest mass 
accretion rate among r <^rB simulations (see Table [2]). 

To assess the robustness of this flaring, we performed 
several checks. Similar flaring is much less pronounced in 
the MHD simulations and in simulations with an isotropic 
thermal conductivity, even though the reduction in M rel- 
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Figure 13. Time and angle averaged density (upper panel) and 
temperature (lower panel) as a function of radius for S3 (oc = 2) 
and MHD3 compared to the initial solutions; the profiles for S3 
and MHD3 are angle and time- averaged. In both S3 & MHD3, the 
plasma is heated to nearly the virial temperature and the density 
in the inner part of the flow decreases significantly. As shown in 
Fig. 1161 there is a significant asymmetry between the polar and 
equatorial regions; the density and temperature are larger at the 
equator while magnetic pressure dominates at the poles. 



Figure 14. The ratio of the time and angle averaged Mach num- 
ber (radial velocity divided by isothermal sound speed Cs) at late 
times as a function of radius for the initial Bondi solution (solid 
line), S2 (dotted line), S3 (short dashed line), and MHD3 (long- 
dashed line). For the calculations with significant magnetic dis- 
sipation (S3 & MHD3), the average infall speed is subsonic at 
r > few rg, and the plasma is in rough hydrostatic equilibrium. 
Although not shown, we find that the radial infall velocity in the 
polar regions is larger than in the hotter equatorial regions. 



ative to Mb is similar; in addition, the flaring is more 
pronounced for larger values of Oc. To assess whether the 
flaring is a numerical artifact, we carried out simulations 
impl ementing anisotropic co nduction using the method of 
iParri sh & Stone (200i) and ISharma et all ([2006), which 
uses simple finite differen cing (our method, based on 
ISharma fc Hammett 1 120071 . limits the heat fluxes so that 
the temperature is always positive). We find quite similar 
flaring with the two methods|j We also changed the magni- 
tude of the artificial viscosity and the explicit resistivity, and 
the form of the explicit resistivity (e.g., 77 ex dr{ci + Va)^ 
instead of eq. |6] where Ca and Va are the local isothermal 
sound speed and Alfven speed, respectively), but the flaring 
in run S3 (Fig. I10|l persisted. These checks argue against 
the flaring being a numerical artifact. The physical origin of 
this flaring is, however, not clear. The MTl could in princi- 
ple be present because the inflow is subsonic, but it should 
be stabilized by magnetic tension. The flaring we flnd is 
qualitatively similar to the "relaxation oscillations" seen in 
some ID time-depend ent models of sphe rical accretion with 
radiative heating (e.g. JCowie et al.ll 19781 ) and could be anal- 
ogous, with spatially extended heating and redistribution of 



^ At late times in the simulations with conduction implemented 
using simple finite differencing, the temperature can become neg- 
ative, which demonstrates the importance of using a method 
th at guarantees that heat flow s from hot to cold, such as that 
of lSha.rma fc Hammetti i2Q0l[) . 



energy by anisotropic conduction precluding the existence 
of a steady state solution. 

Figure [TJ] shows a number of quantities characterizing 
one of the accretion/heating flares over a small interval in 
time. The gravitational potential energy released by accre- 
tion (GMMin/rin) peaks slightly earlier in time than the 
volume integrated resistive heating rate (Hrcs); this is con- 
sistent with the fleld being amplified by fiux freezing and 
then dissipated by resistive diffusion. Figure [T5] also shows 
the maximum values of the local resistive heating rate {r)J^), 
the internal energy per unit volume (e), and the tempera- 
ture, and the minimum value of the density (where the max- 
imum/minimum are over all grid cells in the computational 
domain). The maximum resistive heating rate and maximum 
internal energy have a time-dependence similar to that of 
the accretion energy. The peak temperature and minimum 
density, on the other hand, occur after the majority of the 
resistive heating has taken place, and their evolution in time 
is much more rapid. In fact, the increase in temperature (de- 
crease in density) occurs on a timescale comparable to the 
free-fall timescale at small radii. 

Figure [16] shows contour plots of temperature relative 
to the virial temperature and local resistive heating rate (in 
arbitrary units) in the inner ISr^ at the time of the peak 
in max[T] in Figure [151 (f2ini — 510). Also shown are the 
velocity unit vectors. The resistive heating is strongly con- 
centrated in the equatorial region at < 5rg (near the sonic 
point; see Fig. [14]) where rapid infiow brings together oppo- 
sitely directed field lines. Strong heating leads to modestly 
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Figure 15. Zoom-in on the time-dependence of an accre- 
tion/lieating flare for run S3; note tiiat time liere is in units of 
0.7^ = (GM/r?^)~i/2 not n~^\. The maximum values of the re- 
sistive heating rate {^qj"^), internal energy per unit volume (e), 
and temperature (T), along with the minimum value of the den- 
sity (p), the energy released by accretion (GMMi„/rin), and the 
volume integrated resistive heating rate (Hrcs) are all shown as a 
function of time (the max[X] and min[X] are taken over all grid 
points in the computational domain). All quantities are in arbi- 
trary units. As shown in Figure 1161 the peaks in temperature, 
heating rate, etc. occur at small radii near the equatorial plane. 



Log,„(T/T, 




Figure 16. Left: Login of plasma temperature normalized to 
the virial temperature (T/T^i^) in the inner 15 rg for run S3 
at flint = 510 (the peak in max[T] in Fig. I15I I. Arrows show 
the direction of the velocity vectors. Right: Logio of the local 
resistive heating rate, rjj^ (in arbitrary units), at the same time. 
Both the resistive heating and temperature peak at small radii in 
the equatorial region. 



supervirial temperatures at small radii. This energy diffuses 
to large radii along a finite set of field lines because of the 
anisotropic nature of conduction. The supervirial temper- 
atures also reverse the inflow along particular fleld lines, 
leading to a weak mass outflow (as shown by the velocity 
unit vectors in Fig. I16|l . 



5 DISCUSSION 

We have studied the effects of anisotropic thermal conduc- 
tion on the dynamics of spherical accretion flows using global 
axisymmetric simulations. Our simulations address two dif- 
ferent physical effects that operate on different radial scales 
in spherical accretion flows. 



5.1 The MTI in Spherical Accretion Flows 

We have calculated the nonlinear saturation of the magne- 
tothermal instability (MTI) at radii ~ rs in spherical ac- 
cretion flows (§3), where rs — GM/(?s is the gravitational 
sphere of influence of the central object (of mass M) and Cs 
is the isothermal sound speed of the ambient plasma. The 
MTI does not grow appreciably at radii <^ tb for the canon- 
ical Bondi spherical accretion solution because the inflow is 
supersonic. Our simulations are the first global simulations 
of the MTI, but our results are largely consis t ent with th e 
previous local simulations of IParrish fc Stone I (120051 . 120071 ). 
We have focused on studying the saturation of the MTI for 
initially weak vertical magnetic fields (/3 3> 1). 

We find that the MTI saturates by aligning the mag- 
netic field lines with the temperature gradient, i.e., by 
making the field lines primarily radial (Figs. |2] - |4)| . This 
occurs on approximately the local dynamical time (~ 
[g|dlnr/(iz|]~^'^). Since the field lines end up largely ra- 
dial, the rearrangement of the magnetic field results in a 
radial heat fiux similar to the field-free (Spitzer) value; in 
the saturated state, energy transport by conduction domi- 
nates transport by the convective motions associated with 
the MTI (Fig. [7| . Given the significant radial heat flux, the 
density and temperature profiles of the accretion flow adjust 
to satisfy V • Q ~ 0, where Q is the conductive flux of en- 
ergy. In all cases of Hi{r) that we considered (see H2.1|l . the 
mass accretion rate in the presence of the MTI differs only 
modestly from the canonical Bondi rate (Fig. |6} . 

Our results indicate that the saturation of the MTI is 
quasilinear in the weak fleld limit; the magnetic and con- 
vective energy densities in the saturated state are roughl 
proportional to the initial magnetic energy density (Fig. [S| 
The net result is that the MTI amplifies the magnetic field 
by a factor of ~ 10 — 30 independent of the initial magnetic 
field strength (for weak fields). For stronger initial fields, 
magnetic tension can saturate the MTI before the field has 
become largely radial l|Parrish fc Stone Il2005l . l2007l ). In this 
case, the amplification of the magnetic field will be less than 
we have found here. 

The quasilinear saturation of the MTI can be un der- 
stood qualitatively using linear theory (|Quataert II2008I ); as 



^ Ian Parrish has found this same result in local three- 
dimensional Cartesian simulations (private communication). 
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a result, we believe that this is a general property of the 
MTI, rather than being specific to the spherical accretion 
problem considered in this paper. The linear growth rate of 
the MTI is reduced by a factor of ~ Bg/B for nearly ra- 
dial fields, and the growth occurs on progressively smaller 
scales for small Be/B (eq. [8]). Thus, as the field lines be- 
come radial, the MTI effectively shuts off its own growth. 
This is in stark contrast to the magnetorotational instability 
in differentially rotating plasmas, which is an exact nonlin- 
ear solution of incompressible MHD, and for which growth 
conti nues unimpeded even when 5B ^ B (| Goodman fc Xul 
1 19941 ). 

One simplification in the current calculations is that 
we have not included anisotropic viscosity; the diffusion co- 
efficient for anisotropic viscosity in a collisional plasma is 
smaller than the diffusion coefficient for anisotropic conduc- 
tion by a factor of ~ {mp/nisY''^ ~ 40; in a coUisionless 
plasma, the relative importance of these two effects is more 
complex and depends on small-scale kinetic microinstabili- 
ties (e.g., Sharma et al. 2007). Local simulations of the MRI 
find that the saturation of the instability is sensitive to the 
ratio of the resistivity to the viscosity (e.g., Fromang et al. 
2007). In future calculations, the effects of anisotropic vis- 
cosity on the evolution of the MTI should thus be studied, 
but this is probably first best done in a local calculation, 
rather than in the global geometry considered here. 



5.1.1 Application to the Galactic Center 

One potential application of the MTI at large radii in spher- 
ical accretion flows is to the observed Faraday rotation from 
the Galactic center. Millimeter observations of Sgr A* find 
a rotation measure (RM) of ~ —6 x 10^ rad m~^ that has 
been r oughly constant ove r the past ~ 8 years ([Aitken et al. I 
l200d : iBower et alH l2005l : iMarrone et alTI |2007^ . Propaga- 
tion through the interstellar medium cannot account for this 
large RM. In addition. X-ray observations of the thermal 
plasma in the central parsec of the Galactic center find that 
the electron number density is ~ 100 cm^'^ and the temper- 
ature is « 2 keV at a distance of ~ 10^ rt, f rom the black 
hole l|Baganoff et al. I l2003l : lOuataert I [20ci3 '). If the mag- 
netic field on this scale were uniform, radial, and in rough 
equipartition with the thermal pressure, B ~ 2/3~^" mG 
and RM ~ 5000/3"^'^ rad m~^, two orders of magnitude 
smaller than what is observed. This shows that the observed 
RM must be produced at distances <^ lO^r^ from Sgr A*. 

Our calculations show that given a sufficiently coherent 
seed magnetic field at large radii, the MTI will generate a 
relatively coherent radial field at distances of ~ rs ~ 10^ rg 
from Sgr A*. Given the measured density and temperature 
of the plasma, the electron mean free path is ~ tb at radii 
~ Tb and thus the MTI grows on a dynamical time. Numer- 
ical simulations of accretion onto Sgr A* from stellar winds 
located at ~ 10^ — lO^r^ find tha t the circularization radius 
of the flow is rcirc ~ 10^ - lO^rg l|Cuadra et aLlbOOd l. If the 
Bondi solution is applicable from the scale of the observed 
plasma at ~ rs to ~ rcirc, then dRM/dlnr ex r"^'* at flxed 
/3 and the RM produced by the roughly spherical inflowing 
plasma at large radii around Sgr A* is 



It is difficult to determine with confidence whether the stel- 
lar winds feeding Sgr A* - or the shocks between such 
winds - produce a sufficiently large and coherent 'seed' mag- 
netic field at ~ Tb to account for the RM towards Sgr A*. 
However, for an initial value of /3 of /3o (assumed ^ 1) at 
r_B ~ 10^ rg, the combined action of the MTI and fiux freez- 



-3/2 



at Tcirc . 



RM ~ lO'^ /3 



'1/2 



/ ^circ 



-7/4 



radm 



(13) 



ing will decrease /3 to ~ 10~ /3o('"circ/10 Vg 
Thus even a relatively weak magnetic field (/So ~ 10^ — 10*; 
-Bo ~ 0.1 — l/iG) on scales of ~ lO'^rg ~ 0.1 pc will be am- 
plified by the MTI and fiux freezing to a value capable of 
explaining the observed RM from Sgr A*|j Our explanation 
for the observed RM also naturally accounts for its stability, 
since the dynamical time at ~ 10'^ — 10*rg is ~ 0.3—10 years. 
By contrast, if the RM is produced in the rotating accretion 
flow even closer to the black hole, it is harder to account f or 
the observed stability (ISharma. Quataert, fc Stone 1120071 ). 



5.2 The Effects of Thermal Conduction on the 
Dynamics of Spherical Accretion 

The second set of simulations we have carried out address 
the effects of anisotropic conduction and resistive dissipation 
of magnetic energy on the dynamics of spherical accretion 
flows at small radii ~ 2 — 256 Vg (jj4). Our results are in rea- 
sonable agreement with those of llgumenshchev fc Naravan I 
(2002), who studied the same physical problem but without 
anisotropic conduction. Although the MTI does not dramat- 
ically modify the dynamics of spherical accretion at ~ r_g , it 
does generate a signiflcant radial magnetic field that can be 
further amplified by fiux freezing. For r <C r_g, supersonic 
inflow in (magneto)hydrodynamic Bondi accretion leads to 
magnetic field amplification by flux freezing, converging field 
lines, and resistive heating of the plasma. When the mag- 
netic energy density is comparable to the gravitational en- 
ergy density, resistive heating can heat the plasma up to 
the virial temperature an d significantly modify the dynam - 
ics of spherical accretion (|lgumenshchev fc Naravan II2002I '). 
This is only likely to occur for sub-Eddington accretion rates 
onto compact objects, since only under these conditions is 
radiative cooling negligible. 

All of our calculations with resistivity and strong mag- 
netic fields at ~ 100 rg - both with and without anisotropic 
conduction - show a significant decrease in M after a few 
free-fall times, with M ~ O.OlA/s in the quasi steady state 
at late times (Fig. llOp . The high temperatures generated by 
resistive heating stifie the infiow, reducing the accretion rate 
(Fig. nop and generating a subsonic (Fig. I14p nearly hydro- 
static atmosphere around the central object. For realistic 
values of the thermal conductivity for coUisionless plasmas 
(see >12.ip . we find that conduction transports the majority 
of the energy released by accretion to large radii (Figs. 1111 
fc ll2p . By contrast, although turbulent motions are present, 
there is no net convective energy flux from small to large 
radii, even in our simulations without conduction (Fig. I12prl 
This indicates that signiflcant radial transport of energy (be 



^ If most O stars have surface magnetic fields comparable to 
the ~ kG field m easured in the young O star 9^ Orinis C 
IIDonati et al.ll2002l) , stellar winds alone will fill the central ~ 0.1 
pc with a magnetic field ~ mG. 

^ Our interpretation of the MHD simulations of spheri- 
cal accretion without conduction thus differs from that of 
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it conductive or convective) is not necessary to suppress ac- 
cretion onto the central object; instead, heating that raises 
the temperature to the virial temperature appears sufficient. 
Indeed, it is well-known that there is no supersonic Bondi 
solution for adiabatic indices 7 > 5/3. Heating via resistive 
dissipation and/or artificial/numerical viscosity effectively 
increases the adiabatic index above this value, thus stifling 
accretion onto the central object. Thermal conduction can 
further reduce the accretion rate by conducting the dissi- 
pated energy to larger radii. 

Overall, our resu lt s are in reasonable agreement with 
Ijohnson fc Quataert I l|2007l )'s simplified one-dimensional 
models for the structure of spherical accretion in the pres- 
ence of conduction and plasma heating. Our simulations con- 
firm their hypothesis that a significant fraction of the energy 
released by accretion is transported to large radii by ther- 
mal conduction. In addition, their prediction that M /Mb 
should decrease by ~ 2 — 3 orders of magnitude - with a 
factor of ~ 10 suppression due to heating alone (their Fig. 
8) - is consistent with our numerical results. 

The reduction in M that we find is similar to that 
needed to explain the low lumi nosity of Sgr A* (e.g., 
ISharma. Quataert. fc Stone I |2007|). It is also reasonably 
consistent with t he limits on M from F araday rotation mea- 
surements (e.g., iMarrone et al. Il2007h . However, it is un- 
likely that the present calculations are quantitatively ap- 
plicable to Sgr A* given the estimated circularization ra- 



dius of the accretio n flow in Sgr A* of rdrc 



10'' 



10V„ 



l|Cuadra et al.ll2006l ). The gravitational potential energy re- 
leased accreting to ~ rdrc is unlikely to significantly modify 
the accretion flow dynamics. Thus the accretion rate and 
outward conductive energy flux will be largely determined in 
the rotationally supported flow at radii <_ rdrc Nonetheless, 
our results demonstrate that even in the absence of angular 
momentum, the Bondi accretion solution does not apply to 
dilute magnetized accretion flows. Instead, the accretion fiow 
is heated to nearly the virial temperature by resistive dissi- 
pation of magnetic energy, the accretion rate is <C Mb, and 
most of the gravitational potential energy released by ac- 
cretion is transported to large radii by thermal conduction. 
In addition, in our simulations with anisotropic conduction, 
the accretion and heating are time-dependent, with episodic 
"fiaring" (Figs. fTUlllip followed by mass outfiow along a sub- 
set of the field lines (Fig. I16p . It is important to determine 
if analogous effects are present in rotating accretion flows 
(which is by no means clear), given the possible application 
to the br oad-band electromagn etic flares observed from Sgr 
A* ('e.g.. lMarrone et aLll2008l ) and other systems. 



cooling catastrophically. At larger radii, the temperature 
in clusters increases from r^ ~ 0.1 kpc to ^ 100 kpc, 
where it begins to decrease (|Piffaretti et al.ll2005l ). Thus the 
plasma at intermediate radii in clusters app ears stable to the 



MTI (it is, however unsta ble to the HBI: lOuataert II2OO8I: 



IParrish fc Quataert Il2008h . However, IChandran fc Dennis I 
( 20061 ) have shown that cosmic rays from a central active 
galactic nucleus (AGN) can make the plasma in the cores 
of clusters unstable to a cosmic-ray mediated version of the 
MTI a dpcR/dr + nkdT/dr < {pcR is the cosmic ray pres- 
sure), i.e., if there is a sufficiently large cosmic ray flux from 
the central AGN. It is likely that the MTI in the presence of 
cosmic rays saturates in a manner analogous to that found 
here, by aligning the field lines primarily in the radial di- 
rection, with a roughly constant radial cosmic ray and con- 
ductive luminosity. One of the outcomes of this work (see 
also Parrish & Stone 2007!) is that convective transport of 
energy is subdominant in the presence of the MTI; this may 
rule out turbulent heating due t o cosmic-ray driven convec- 
tion (jChandran fc Rasera Il2007h as an important source of 
plasma heating in cluster cores. However, if the cosmic-ray 
mediated MTI generates a largely radial magnetic field in 
cluster cores, the conductive energy fiux from large radii 
can produce significant heating at small radii, analogous to 
what we have found in our global simulations (see Fig. [?])• 
It is clear that a careful study of the combined action of the 
MTI and HBI in the presence of cosmic-rays is necessary to 
understand the thermal structure of galaxy clusters. 
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5.3 Additional Applications 

In addition to the applications considered here, thermal con- 
duction is also important in the intracluster medium of 
galaxy clusters IjSarazin 1 1 19861 ). A conductive energy fiux 
from the accretion fiow onto the central black hole in the 
nucleus of galaxy clusters could help heat the interclus- 
ter plasma at radii ^ 0.1—1 kpc and prevent it from 



llgumenshchev fc Naravan I l l2002h . who emphasized the impor- 
tance of radial energy transport by convection. This needs to be 
investigated further. 
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